def llt_pred_interval(y_pred, t, F, f, var, alpha=0.05):
return y_pred + np.array([1, -1]) * stats.t.ppf(
1 - alpha / 2, t - param_n
) * np.sqrt(var) * np.sqrt(1 + f.T @ np.linalg.inv(F) @ f)
def f(t):
return np.array(
[
[
1,
t,
np.sin(2 * np.pi / p * t),
np.cos(2 * np.pi / p * t),
]
]
).T
def llt_predict(start=10, lamb=0.9):
assert start > 0, "Start must be greater than 0"
# Createnp.linalg.inverse transition matrix
L_inv = np.linalg.inv(
np.array(
[
[1, 0, 0, 0],
[1, 1, 0, 0],
[0, 0, np.cos(2 * np.pi / p), np.sin(2 * np.pi / p)],
[0, 0, -np.sin(2 * np.pi / p), np.cos(2 * np.pi / p)],
]
)
)
# Create design matrix for t'th time step
f_0 = f(0)
F = np.zeros((param_n, param_n))
h = np.zeros((param_n, 1))
# Store predictions and variances
y_pred_train = []
y_pred_train_interval = []
y_pred_train_mean = []
y_pred_train_var = []
y_pred_train_var_sum = 0
thetas = []
# Go through each data point and make one-step predictions
for i in range(0, train_n):
# If burn-in period is over, make predictions
if i >= start:
F_inv = np.linalg.inv(F)
# Calculate parameters for previous time step
theta = F_inv @ h
# Make one step prediction using parameters
# from previous data to predict current data.
y_pred = f(1).T @ theta
# Calculate noise variance from next prediction
y_pred_train_var_sum += (
(y_train[i, 0] - y_pred[0]) ** 2 / (1 + f(1).T @ F_inv @ f(1))
)[0, 0]
y_pred_train_var.append(y_pred_train_var_sum / (i - start + 1))
# Store predictions
y_pred_train.append(y_pred[0, 0])
y_pred_train_mean.append(theta[0, 0])
thetas.append(theta)
# Store prediction interval
y_pred_train_interval.append(
llt_pred_interval(y_pred, i - 1, F, f(1), y_pred_train_var[-1], 0.05)
)
# Update F and h
F = F + lamb**i * f(-i) @ f(-i).T
h = lamb * L_inv @ h + f_0 @ y_train[i : i + 1]
# Calculate final parameters
theta = np.linalg.inv(F) @ h
thetas.append(theta)
# Make predictions for test data
y_pred_test = np.array([f(t + 1).squeeze() for t in range(test_n)]) @ theta
# Make intervals for test data
y_pred_test_interval = np.array(
[
llt_pred_interval(
y_pred_test[t], train_n, F, f(t + 1), y_pred_train_var[-1], 0.05
)
for t in range(test_n)
]
).squeeze()
# Return a big mess
return (
np.array(y_pred_train),
np.array(y_pred_train_interval).squeeze(),
np.array(y_pred_train_mean),
np.array(y_pred_train_var),
y_pred_test,
y_pred_test_interval,
)
# Amount of data to skip predictions for
burn_in_1 = 10
# Retrieve predictions
(
llt_y_pred_train,
llt_y_pred_train_interval,
llt_y_pred_train_mean,
llt_y_pred_train_var,
llt_y_pred_test,
llt_y_pred_test_interval,
) = llt_predict(start=burn_in_1, lamb=0.9)
